function [forc_SOn,forc_Soff,forcDecomp]=zb_quickforecast_shockdecomp(forc_pos,GG,RR,ZZ,CC,nforc,sT,etaMat,pos_gdpdef,in,shonames,gstruct,groupnames) 
% =========================================================================
% =========================================================================
% ZB_QUICKFORECASTS_SHOCKDECOMP
% 
% Generate model-based forecasts. Same as ZB_QUICKFORECASTS but also
% decompose across the different shocks in etaMAT which is the matrix of
% shocks that must have as many rows as NFORC
% 
% function [forc_mat,yforc]=zb_quickforecast(forc_pos,GG,ZZ,CC,nforc,sT,pos_gdpdef,popsTF,in) 
% 
%% 1. Inputs 
% -------------------------------------------------------------------------
% GG,ZZ,CC  :Solution matrices of model 
% NFORC     :Number of Forecasts    
% ST        :Last data state, e.g. ST|T
% etaMAT    :Matrix of future shocks, e.g. ETA(T+1),
% ETA(T+2),.....,ETA(T+nforc). Matrix will *NOT* be padded inside the code 
% fill with zeros outside if needed 
%
% POS_GDPDEF:Position of GDPDEF *in STANAMES* 
% POPST     :Last data of population growth (annualized) matching end of sample 
% POPF      :Series of forecast population growth, matching forecast period 
%            *NOTE* if length(POPF) is less than NFORC, last point will be
%            extended, but user will be queried about it. 
%
% IN requires fields
% .CONSNONADJ:      Non-adjusted constant 
% .ADDP:            Prices flag 
% .POP :            Population flag 
% .PARAM_GDPSS:     Constant for GDP Deflator 
%
%% 2. Outputs 
% -------------------------------------------------------------------------
% *Note:* last state is in the first row 
% FORC_SON:     [length(forc_pos) NFORC] Matrix of forecasts with all
%               shocks on 
% FORC_SOFF:    [length(forc_pos) NFORC] Matrix of forecasts with all
%               shocks off
% FORC_DECOMP   [length(forc_pos) NFORC NGROUPS] Decomposes the difference
%               between forecasts on and off across groups, and check these
%               sum. 
% Alejandro Justiniano December 1st 2010 
% =========================================================================
% =========================================================================
consnoadj       =in.consnoadj ;
forc_scale      =in.forc_scale;
forc_addp       =in.forc_addp;
forc_pop        =in.forc_pop ;
param_gdpss     =in.param_gdpss;
in.forc_pos=forc_pos;
in.pos_gdpdef=pos_gdpdef;


%% Padd population series with last entry if necessary 
% if length(popsTF) < nforc
%     disp('Extending Population series with end-point to cover forecast horizon'); 
%     quer('c'); 
%     paddf=popsTF(end)*zeros( nforc - length(popsTF) ,1  ); 
%     popsTF =[popsTF;paddf]; 
% end 

%% Check ETAMAT is of right size 
if size(etaMat,1)~=nforc 
    error('ETAMAT must have as many rows as forecasts') 
end 

%% forc_SOFF: Generate Forecasts with ALL SHOCKS set to zero 
% This is the unconiditional forecast. 
% sT is unstransformed final state, in deviations from model mean 
% Number of actual forecasts is NFORC; last state is NOT first entry 
nser=length(forc_pos); 
[stforc_Soff]=zb_forecast_loop(GG,RR,ZZ,CC,sT,[],nforc); 
forc_Soff=stforc_Soff(:,forc_pos);
for ii=1:nser;  
    if forc_addp(ii)    ==1;
        forc_Soff(:,ii)  =forc_Soff(:,ii)-stforc_Soff(:,pos_gdpdef);
    end
end


% Following lines replaced by function ZB_FORECAST_LOOP 
%  checked identical 11/18/11
% stforc=zeros(nst,nforc); 
% yforc =zeros(size(ZZ,1),nforc); 
% stforc(:,1)=GG*sT(:); 
% yforc( :,1)=ZZ*(CC+stforc(:,1) ); 
% for ii=2:nforc; 
%     stforc(:,ii)=GG*stforc(:,ii-1);
%     yforc(:, ii)=ZZ*(CC+stforc(:,ii) ); 
% end
% stforc=stforc'; 
% yforc =yforc' ; 

%% Tranform forecast 

%% Following lines replaced by function ZB_QUICKFORECAST_SUB.m 
%  checked identical 11/18/11
% Obtain Annualized GDP 
%forc_SOff_mat=zb_quickforecast_sub(stforc_Soff,popsTF,in); 
% if any(forc_addp ~= 0 ); 
%         gdpdef_fancons=4*( stforc(:,pos_gdpdef)+param_gdpss ); 
% end 
% % Transform forecasts 
% nser=length(forc_pos);
% forc_mat=zeros(nforc,nser);
% for ii=1:nser;
%     forc_mat(:,ii)      =forc_scale(ii)*(stforc(:,forc_pos(ii))+consnoadj(ii) ) ;
%     if forc_addp(ii)    ==1;
%         forc_mat(:,ii)  =forc_mat(:,ii)-gdpdef_fancons;
%     end
%     if forc_pop(ii)==1;
%         forc_mat(:,ii)  =forc_mat(:,ii)+popsTF;
%     end
% end
% % Check subroutine that will be used to generate subsequent forecasts works well 
% forcmat2=forc_mat; 
% forc_mat=zb_quickforecast_sub(stforc,popsTF,in); 
% comparemat(forc_mat,forcmat2); 


%% stforc_S0n: Generate forecast with all shocks ON 
stforc_SOn=zb_forecast_loop(GG,RR,ZZ,CC,sT,etaMat,nforc); 
forc_SOn=stforc_SOn(:,forc_pos);
for ii=1:nser;  
    if forc_addp(ii)    ==1;
        forc_SOn(:,ii)  =forc_SOn(:,ii)-stforc_SOn(:,pos_gdpdef);
    end
end


%% Extract groups 
gsum=zb_extractgroups(gstruct,shonames,groupnames); 
ngroups=length(groupnames); 
forcDecomp=zeros(nforc,length(forc_pos),ngroups); 
for ii=1:ngroups 
    pos=getfield(gsum,['p',num2str(ii)]); 
    eta_temp=zeros(size(etaMat) ); 
    eta_temp(:,pos)=etaMat(:,pos); 
    st_temp=zb_forecast_loop(GG,RR,ZZ,CC,sT,eta_temp,nforc); 
    forc_temp=st_temp(:,forc_pos); 
    for jj=1:nser;
        if forc_addp(jj)    ==1;
            forc_temp(:,jj)  =forc_temp(:,jj)-st_temp(:,pos_gdpdef);
        end
    end
    forcDecomp(:,:,ii)=forc_temp-forc_Soff; 
end 
maxdif=comparemat(squeeze(sum(forcDecomp,3)),forc_SOn-forc_Soff);
if maxdif > 1e-4;
    disp('Conditional Forecasts cannot be decomposed across shock groups');
    error('Conditional Forecasts cannot be decomposed across shock groups');
end 
